python实现梯度下降->牛顿->拟牛顿算法

本文主要抄自这篇文章, 我只是学习一遍, 然后加入了自己的理解, 并不是比他写的更好. 因为我是做开发的, 所以文章会更偏工程.

首先澄清, 我们要解决的问题是最小化下面这个方程:

\[ \min_{x\in\mathbb{R}^2}f(x)=100(x_1^2-x_2)^2+(x_1-1)^2 \]

思路

梯度下降法是求函数最小值的常用方法, 思路是, 求函数f(x)的梯度g(x), 梯度的作用是, 可以指明求解x的方向, 因为梯度的反方向就是函数值下降的方向. 所以, 我们任意设置一个初始值\(x_0\), 然后让\(x_0\)做如下迭代:

\[ x_1 = \alpha g(x_0) + x_0 x_2 = \alpha g(x_1) + x_1 \cdots x_(k+1) = \alpha g(x_k) + x_k \]

经过这个过程, 函数f(x)的值会逐渐降低. 这个过程有两个未知的变量\(\alpha\)和\(g(x)\), 前者被称为步长, 它越大, 迭代过程x的变化越快, 但是可能步子太大跨过最优点; 后者是梯度, 也就是自变量的偏导数构成的向量. 一步一步来, 不用多想, 先求得f(x)的梯度g(x).

梯度

函数f(x)的梯度如下所示:

\[ g(x)=\begin{pmatrix} 400x_1(x_1^2-x_2)+2(x_1-1) \\ -200(x_1^2-x_2) \end{pmatrix} \]

但是, 我很同情你们不懂微积分的程序员们, 我给你指条明路, 使用python的sympy模块可以进行符号运算, 对于一些简单的函数可以快速得到导数/偏导数, 比如我们要求上面目标函数f(x)的梯度.

求f(x)关于\(x_1\)的偏导数:

1
2
3
4
5
6
7
8

from sympy import *
x1 = Symbol('x1')
x2= Symbol('x2')
f=100*(x1**2-x2)**2+(x1-1)**2
# 关于x1 求偏导数
diff(f, x1)
# 关于x2 求偏导
diff(f, x2)

输出:

\[ \begin{equation*}400 x_{1} \left(x_{1}^{2} - x_{2}\right) + 2 x_{1} - 2\end{equation*} \]

\[ \begin{equation*}- 200 x_{1}^{2} + 200 x_{2}\end{equation*} \]

不过, 对于比较复杂的函数, sympy可能会卡死, 所以不要指望sympy能帮你跨过数学这道门槛.

求α

我们可以很任性, 设置α为一个固定常量, 比如 0.1. 后面再讲选择α的高级方法.

设置任意的初始\(x_1, x_2\)

就是这么任性, 将两个自变量的初始值设置为(7, 7)

迭代

根据上面的思路, 我写了下面的python代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

# 目标函数
f = lambda x:100*(x[0]**2-x[1])**2 + (x[0]-1)**2
# 梯度函数
g = lambda x:([400*x[0]*(x[0]**2-x[1])+2*(x[0]-1), -200*(x[0]**2-x[1])])
# 初始x0
x0 = [7, 7]
# 最大迭代次数, 防止不收敛, 进入死循环
max_loop = 5000
# 迭代终止条件
epsilon = 1e-5
# 步长
alpha = (0.1, 0.1)
# 迭代计数器
k = 0
# 第一个自变量x1 的偏导数

# 迭代的x
x = x0[:]

while k < max_loop:
    delta = g(x)
    value = f(x)
    print(delta)
    for i in range(2):
        x[i]=alpha[i]*delta[i] + x[i]
    print('new x:', x)
    if abs(f(x)-value) < epsilon:
        break

print('Loop:', k)
print('x:', x)

参考资料

• http://blog.csdn.net/golden1314521/article/details/46225289